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Abstract 



The MACRO underground detector at Gran Sasso Laboratory 
recorded 60 million secondary cosmic ray muons from February 1989 
until December 2000. Different techniques were used to analyze this 
sample in search for density excesses from astrophysical point-like 
sources. No evidence for DC excesses for any source in an all-sky 
survey is reported. In addition, searches for muon excess correlated 
with the known binary periods of Cygnus X-3 and Hercules X-1, 
and searches for statistically significant bursting episodes from known 
gamma ray sources are also proved negative. 



1 Introduction 

The interest in the search for cosmic ray point sources identified from a 
measurable flux of underground muons has historical motivations mainly 
because of the "CygX-3 saga". Cygnus X-3 (CygX-3) is a galactic binary 
system well studied in all types of electromagnetic radiation, most notably 
in the X-rays. At 7-ray energies, E > 10^^ eV, the Kiel Extensive Air 
Shower array (EAS) initially reported an excess of events from CygX-3 
correlated with its 4.8 hour binary period [Q. This observation appeared to 
confirm previous results at lower energies 0, |^ . Subsequently, several groups 
operating EAS array experiments and atmospheric Cerenkov telescopes 
confirmed the signals with different statistical significance and at different 
energy thresholds [^, ^, H, ^ ||, . These results prompted a new generation 
of more sensitive experiments using new techniques. The new experiments 
reported clear evidence of TeV 7-ray emissions from many galactic {i.e., Crab 
and SN1006 [U], |12[) and the extragalactic AGN-like (i.e., Mrk421 and 



MrkSOl [|T^, [1^) sources. Finally, in 1990-1991 the new generation detectors 
CASA-MIA[ll5l and CYGNUS put stringent upper limits to the 7-ray 
signal from CygX-3, which excluded early observations. 

Soon after the initial EAS detection from CygX-3, two underground 
experiments (Soudan 1 |]T^ and NUSEX []TB[) reported excesses in TeV muons 



pointing back to the CygX-3 direction when correlating event arrival time 
with the known orbital period of the source. These detections suggested new 
physics beyond the Standard Model since muons are the products of the decay 
of charged pions and kaons, which are produced by primary Cosmic Ray (CR) 
interactions with atmospheric nuclei. Charged CR nuclei cannot propagate 



directly from CygX-3 across 10 kpc to the Earth and there is only a very 



low probability of UHE 7-rays producing TeV muons underground M, 20 



Investigations by other underground detectors pT], |23, 53, E3, however failed 



to detect any significant muon excess from that source or any other. 

The MACRO experiment, which ended its operations life in December 
2000, was a large area underground detector able to reconstruct muon arrival 
directions to very high accuracy The apparatus started operations in 
1989 and was completed in April 1994. A preliminary investigation using a 
limited data sample of ~ 1.8 x 10^ muons collected by its first and second 



supermodules was published in 1993 p6[. Since then, MACRO has increased 



the muon sample by a factor of 30, its pointing capability has been accurately 
determined and its direction reconstruction capability carefully studied. 

In this paper we present the final results on the search for a muon excess 
with respect to the evaluated background by surveying the sky in declination 
from -15° to 90° and from several candidate UHE gamma ray emitters. 
We also present a more sensitive search for CygX-3 and HerX-1 using their 
known periodicities. Finally we searched for bursting behavior of Mrk421 
and MrkSOl. 



2 The MACRO detector 

MACRO was a multipurpose underground detector located in the Gran Sasso 
underground Laboratory (LNGS-INFN). It was designed to search for rare 
events in the cosmic radiation and its sensitivity was optimized to detect 



supermassive magnetic monopoles. The detector ||2^, ^ had a modular 
structure and dimensions of 76.5 x 12 x 9.3 m^ with a total acceptance for an 
isotropic flux of about 10^ m^sr. The rock overburden had a minimal depth 
of 3150 m.w.e and an average of 3700 m.w.e. 

The minimal energy of surface muons that trigger the MACRO apparatus 
was about 1.3 TeV. The rock coverage was very irregular and the actual 
slant depth under which a source is observable was taken into account when 
evaluating the muon flux at the surface. The full detector acceptance for 
downgoing muons with zenith angles < 72° was about 3100 m^sr. 

The detector has worked with different configurations starting with 
1 supermodule at the beginning and finally with the full configuration 



consisting of 6 supermodules and the upper half called the "Attico" p8[ . The 



active detection elements are planes of streamer tubes for tracking and liquid 



scintillation counters for fast timing. The lower half of the detector is filled 
with streamer tube planes alternating with trays of crushed rock. The upper 
part is hollow and contains the electronics racks and work areas. There are 
10 horizontal streamer tube planes in the bottom half of the detector, and 4 
planes on the top, all with wires and 27° strips, providing stereo readout of 
the detector hits. Six vertical planes of streamer tubes cover each side of the 
detector. 

3 Data selection 

The data used in this analysis were collected in the period May 1989- 
December 2000. The total number of recorded muons is approximately 60 
million. 

In order to optimize the quality of the tracking reconstruction and to 
eliminate periods when the detector was malfunctioning, we applied several 
cuts on a run-by-run basis. These run cuts include a check on the efficiency of 
the streamer tubes system whose data were used for the track reconstruction, 
and a check on the event rate. The streamer tube efficiency was obtained 
using the sub-sample of tracks crossing all the 10 lower streamer tube planes; 
runs having an average efficiency smaller than 90% for the wire view and 85% 
for the strip view were discarded. As the detector configuration changed 
during the data taking, the average value of the counting rate was computed 
for each run. Those with abnormal rates, i.e., runs having muon counting 
rate that deviated by more than ±3cr from the average, were cut. 

In addition to the described run cuts, we applied event cuts by excluding: 

• events having zenith angles larger than 72°; the cut is due to the large 
uncertainties in the rock depth crossed by these events; 

• events whose reconstruction in one of the two streamer tube stereo 
views is missing; 

• muons in a bundle with multiplicity larger than 2; this cut is necessary 
because the reconstructed direction in space for high multiplicity events 
is not completely reliable; 

• events that had no arrival time given by the Universal Time Clock 
(UTC); 



events that do not cross at least 3 planes of the lower part of the 



apparatus; this cut excludes low energy muons p9[ coming from 
possible secondary interactions and it improves the angular resolution. 
The minimal muon energy to satisfy this cut is 1 GeV. 

These cuts reduce the sample to 49.9 million well-reconstructed single and 
double muon events during 74,073 hours of livetime. Table |I| shows, for 
different data taking periods, the detector configuration, the number of events 
which survived the cuts and the effective livetime. 



4 Background estimation 

Several methods were suggested to evaluate the background when looking 
for an excess of counts from a fixed direction in the sky ]TB| , Many 
experiments simply consider, as an estimation of the background, the average 
value of the counts of the surrounding sky bins with respect to a selected one 



| 32| . An alternative way is to average the counts of all the bins at the 
same declination 6, except the chosen one in a particular value of the right 
ascension (RA). The assumption that cosmic rays arrive uniformly from any 
direction in the sky is implicit in this method. For an underground detector 
both methods require that the unequal distribution of the overburden and 



the dead-time of the apparatus be taken carefully into account |2T[ . 

In our analysis, we adopted a different approach. Assuming that the 
arrival directions and the arrival time of underground muons are uncorrelated 
variables, we constructed sets of simulated events by randomly coupling the 
times and the directions of each event in local coordinates. A mandatory 
requirement for this method is a good accuracy in the measurement of the 
arrival direction and of the recorded time. For this reason, we excluded from 
the analysis all the runs and the events with errors in the readout time, as 
described in the previous section. We determined that the optimum number 
of background events to be associated with each real event, minimizing the 
computer processing time, is 25. A total number of 1.3 x 10^ simulated muons 
were generated as background. Figure |l] presents the distribution of observed 
events vs. the right ascension RA (a) and the sine of the declination sin5 (b) 
together with the simulated background. The fluctuations in the simulated 
distributions are small compared with those in the data. The unevenness in 
RA is due to the effect of the dead times. The sm6 distribution reflects the 
shape of the mountain and also the exposure of MACRO (see |p6|). 



5 Best sky bin definition 



It is usual for astronomical telescopes to define a point spread function (PSF) 
that in most cases is a simple bi-dimensional normal distribution. The PSF 
can be used to define the optimal source bin, i.e., the half-angle 6* of a cone 
centered on a source giving the maximum signal over background ratio 
For a normal PSF distribution with variance a it is shown in [|l^ that the 
best value of the half-angle of a cone is 6' = 1.58cr. This cone contains, on 
average, 72% of the total number of events. 

The MACRO PSF was defined using the double muon sample. A 
double muon event (as in general multiple muon events) is produced by an 
interaction of a primary nucleus at the top of the atmosphere; muons are 
expected to arrive at the earth surface practically parallel. When reaching 
an underground detector, the reconstructed spatial directions of these two 
muons differ as a consequence of the independent scattering of each muon. 
We can therefore estimate the MACRO PSF by using as a variable the 
distribution of the angular separation of the reconstructed muon directions 
6 for double muon events. This value must be divided by a \/2 factor to 
account for the two independent scatterings . 

We found ||26|, |27|, B8|l that 50% of double muons events are contained 



in a cone of 0.5° half-angle, and that 72% of events are contained in a 
cone of 1.05° half-angle. Initially, this result suggests that the best value 
for the MACRO search bin (assuming that the 9 is normally distributed) 
is an half-angle cone of 1.05°/\/2- To produce the bidimensional {RA,5) 
distribution representing our PSF (shown in Fig. H) we used the differences 
between the reconstructed muon directions in RA and declination 5, divided 
by V2. The main contribution to the differences between the reconstructed 
directions is the multiple Coulomb scattering (MCS) of muons in the rock. 
The PSF shown in Fig. ^ has a non-gaussian shape, with more events for 
small displacements and longer tails. The latter were produced by stochastic 
interaction process in the muon energy loss. 

To obtain the value of the best search cone angle, we performed a 
dedicated Monte Carlo. A source was simulated in a particular {RA, 5) 
cell, and the events were extracted according to the PSF of Fig. ^. The 
background events were generated according to a fiat density distribution. 
We calculated the ratio as a function of the search cone angle, as shown 
in Fig. ^. The maximum of the ratio is reached for a search cone angle 



of ~ 0.4-0.5°, which thus will represent the best choice value for MACRO. 

As an independent check of the validity of the simulation, we studied 
the signal due to the Moon shadowing effect (in this lack of events) 

measured by our detector [^. The cumulative number of missing events 
obscured by the Moon's disk, as function of the angular distance from the 
Moon center, shows maximum evidence for a deficit using a cone of 0.45° half- 
angle. 

From all these indications the best value for the search cone angle was 
chosen to 0.5°. However, in the analyses with very low bin contents such 
as those used in the search for flaring activities, we use an enlarged search 
window due to poor statistics. 



6 Search for DC sources 

In the search for a steady excess of muons from any direction of the sky (DC 
point-like sources search) we performed an all-sky survey without a priori 
assumptions. We divided the sky into 37176 bins of equal solid angle (Af2= 
2.3 X lO^^sr; ARA = 1°, Asin5 = 0.013). These bins have the same AO 
as a narrow cone of half-angle 0.5°. We examined the sky bins looking for 
significant deviations from the simulated background. The deviation was 
defined as a = i2£M^^£££l^ where noi,s is the number of events observed in 
each bin and Uexp is the expected number of background events in that bin. 

We used three different sky grids, each displaced by half bin width in 
RA, in sin((J) and in both coordinates, to take into account the possibility 
that a source can be located at the edge of one of the bins. Fig. ^ shows 
the distributions of the deviations for the four sky-maps with the best-fit 
Gaussian function superposed. The positive deviations in the first map are 
reported in the bidimensional map in RA and sm6 of Fig. ^. No particular 
pattern or clustering of positive deviations is observed. The line indicating 
the galactic plane is also shown. We conclude that there is no evidence for a 
steady source emitting muons in our data. 

We calculated the 95% C.L muon flux for all sky bins using the formula 



Jfn95%)<^^^cm-V^ (1) 



where: 



- n^{95%) is the upper limit for the number of muons in the bin at 
95%C.L., where Uobs and n^xp are respectively the number of observed 
and expected events in that bin. ?7,^(95%) was calculated, according to 
|]33[| , as the value for which: 

2 roo ^-{n^-nj:)'^/2a^ 

7= dn. = 0.05 (2) 

VTT in^{95%) V2a- 

with ni;=nobs - Uexp and a'^=nexp- 

- Aeffii) is the average effective detector area for every bin. It was 
computed by averaging the projected area seen by each muon and 
taking into account the geometrical and the tracking reconstruction 
efficiencies: 

^^ff(') = zrrr^ ^ MRA,s.) (3) 

nobs{i) j=i 

where RAi and Si are the muon arrival right ascension and declination. 

- Texp is the exposure time computed as the time in which each bin is 
visible. A bin is not visible in our apparatus, if its zenith angle is larger 
than 72°. 

- K is a correction factor which takes into account the fraction of muons 
within the bin dimensions for an hypothetical source placed at the bin 
center. We calculated K=0.5 for a search cone of 0.5° half-angle and 
K=0.72 for 1.0° half-angle. 

The distribution of the 95% C.L. upper limit J^*'^2^(95%) for all the 
37176 bins in the sky region accessible for MACRO has an average value 
of 2.3 X 10~^^ cm~^s~'^. For almost all bins, the upper limit ranges between 
1 X 10-13 < j^i'^?/(95%) < 4 X 10-13 cm-V^ 

Finally we also investigated some selected interesting point sources (in 



Table |^) identified by surface telescopes and EAS arrays ^ |3l 

None of the selected sources in the list exhibits a significant deviation from 
the background. 



7 Searches for modulated signals from CygX- 
3 and HerX-1 



Periodic UHE gamma ray sources are attractive from the observational point 
of view. When it is known that the emission from the source is modulated 



with a certain periodicity, the signal-to-noise ratio is improved by yA^ p8 
by folding the event arrival time modulo the source period into N bins. This 
folding procedure was employed to analyze the data from the direction of 
CygX-3 and HerX-1 [|, |9], . 

We used the quadratic ephemeris reported in Table ^ and obtained from 
the fit to the observed X-ray light curve The phase diagrams for 



events coming from a cone of 1°, centered on the position of the two sources, 
are shown in Fig. p. The expected number of background muons in each 
phase bin is also shown with the dashed lines. 

The largest positive deviation (1.8 a) in the CygX-3 cycle is in the phase 
bin 0.1-0.2. The largest positive deviation for HerX-1 (1.7 a) is in the 
phase bin 0.7-0.8. Using the values of the largest fluctuation we computed 
the 95% C.L. upper limits to the modulated muon flux using equation |l|, 
and the results are reported in Table ^. In Fig. |^ our computed upper 
limit for the modulated emission from CygX-3 is compared with those of 
previous underground experiments. Our level is the lowest reached by such 
detectors. The problem of understanding past positive observations still 
remain unsolved. 



8 Search for flaring activity from Mkn421 
and Mkn501 

During the 1995-1999, several Cerenkov telescopes reported observations of 
flaring activities up to 10-20 TeV from the two celestial objects Mkn421 



and Mkn50lQ, These observations prompted us to search for 

bursting muon signals in a 1° half-cone, around the position of these sources 
of U.H.E. photons. We used two different methods. 

In the first method we studied, similarly to |^ , the accumulation rate of 
events coming from each of the two sources, adding day-by-day the differences 
between the measured number of events and the calculated background. 
Figures | and |^ show the cumulative excesses as a function of date (Julian 



days) since MACRO starting data taking. The cumulative excess presents 
fluctuations, but never becomes significant. 

In the second method we assumed (as in [^) that the background has 
a Poissonian distribution. If Uexp is the expected daily background from the 
direction of a source, then the probability to observe a random fluctuation 
of the background large as the observed Uobs events in a day, is given by: 

where a is the ratio of the ON-source time to the OFF-source time. Because 
we extract 25 simulated background events for each real one, a is 0.04. 
For a Poissonian background, the cumulative frequency distribution of P 
is expected to be a power law with index -1, in logarithmic scale. 

Figures 110 and 111] show the value of -logF evaluated day-by-day for 



Mkn421 and MknSOl, respectively. In Fig. |10| the date of the largest 
fluctuations for Mkn421 with respect to the background are also indicated. 

To verify the probability of such a positive fluctuation from Mkn421, we 
computed the quantity -logF for a set of selected bins having an exposure 
similar to the two Mkn objects, each monitored for about 3600 days. 



Figure |T2| shows the log-log plot of the cumulative frequency distribution 
for all these sky bins. Since this cumulative distribution has a slope close to 
-1, as expected in the case of no source detected, we can use the Poissonian 
statistics to compute the expectation to see large fluctuations. 

We set as an attention level a probability of < 10^'^ for a fluctuation. 
We found 4 days with probability larger than 10~^ for Mkn421, and none 
for MknSOl. We observed the Mkn421 source for 3600 days: assuming as 
attention level for any source the probability value of 10~^, the expected 
number of random fluctuations with probability lower than 10^^ is 3600 • 
(1 X 10~^) = 3.6. Therefore the Poissonian probability to observe 4 random 
fluctuation in the same period, with an average of 3.6, is about 20%. This 
probability value therefore excluded a positive observation of a burst from 
this source. 



9 Conclusions 



Since February 1989, the MACRO detector collected 49.9 million well- 
constructed muons. Using this sample, we searched for muon excesses above 



background from all visible sky directions and from known astrophysical 
sources. No significant excesses were found from the all-sky survey. We 
computed the 95% confidence level upper limit J^*'^^(95%) for a steady muon 
flux for all the 37176 sky cells; the average value of J^*'^^(95%) is equal to 
2.3x10"^^ cm~^ s~^. We analyzed the muons coming from the direction of 
Cyg X-3 and Her X-1, searching for a modulated emission, with a negative 
result. The search for a steady or episodic emission coming from Mkn421 and 
MknSOl was made with two different methods. We found no muon excess 
above the estimated background. 
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Period 


Number of 


Number of 


Exposure 




SuperModules 


events 


Time (hr) 


Feb,27 1989-May,20 1989 


1 


244333 


1942.9 


Nov,ll 1989-May,10 1990 


1 


365148 


3072.1 


May,10 1990-Jul,5 1991 


1,2 


1308311 


5274.9 


Jul,5 1991-Apr,29 1994 


1-6 


11549606 


16247.8 


Apr, 29 1994-Dec,15 2000 


1-6+A 


36490390 


47535.4 



Table 1: Description of the muon data sets collected by MACRO. As 
construction of MACRO progressed, more supermodules came on-line and 
increased MACRO'S collection area. Most of MACRO'S muon data was 
collected after MACRO was fully operational in 1994. 
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Tabic 2: Search for muon excesses from selected sources, using half-cones 
of 0.5° and 1°. The number of muons observed Uobs and expected rtexp are 
tabulated, with the quantity a = (2£££z2£sl to indicate the significance of 
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Cyg X-3 |411 


Her X-1 |43|] 


Ephemeris 


P[d] 


0.19968271 
(±2.4x10-"^) 


1.700167788 
(±1.1x10-^) 


P/P [d/yr] 


(1.17 ±0.44) X 10"*^ 


< 2 X 10"« 


To [JD] 


2440949.8989 
(±0.0012) 


2442859.726688 
(±7x10-^) 


Flux 
limits 1° 


jmod (95%) 


1.4 X 10"^^ 


1.6 X 10"^^ 



Table 3: Ephemeris parameters for Cyg X-3 and Her X-1 and modulated 
flux limits computed using equation |I] in a 1° cone around the source position. 
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Figure 1: (a) Right Ascension distribution of the sample of 49.9 miUion 
muons. (b) Distribution of the sin 5 for the same sample of events. The 
normalized simulated background is superimposed in both figures, however 
in (b) they are too close to distinguish. 




Figure 2: MACRO Point Spread Function, derived from the measured 
differences in right ascension and dechnation coordinates of each muon in 
double muon events divided by \/2. 




Figure 3: Simulated function for the MACRO PSF vs the space angle 
from the source center. The maximum value occurs about at 0.45°. 
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Map-2 (a) 
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Figure 4: Map-1 shows distributions of the quantity a = i^E^^-^ssi^ where 
Hobs is the number of events observed in each sky bin and Uexp the expected 
number of background events in that bin. Each bin has the same sohd angle 
An = 2.3 X 10~Sr. Map-2 through Map-4 were obtained by the same 
procedure of Map-1, but with shifts in RA by -|-0.5°(Map-2), or in sin 5 by 
+0.013 (Map-3) or both (Map-4). 
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Figure 5: Bidimensional distribution (in right ascension and sine of the 
dechnation) of the positive value of o" = , where riobs is the number 

of events observed in each bin and n^xp the expected number of background 
events in that bin. The hne superimposed indicates the galactic plane. No 
point sources are seen. 
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Figure 6: Phase diagrams for muon events from 1° half-angle cone centered 
on (a) Cyg X-3 and (b) Her X-1. The dashed histograms represent the 
simulations of the background. 
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Figure 7: Searches for modulated signals from Cyg X-3: present situation for 
the flux limits reported by other experiments |17, |l8[ 21, 22, 23, ^ at 95% 
CL. For MACRO, the method from ref.|3^ was used. 
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Figure 8: Cumulative muon excesses from the direction of Mkn421 (1 "half- 
angle) . 
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Figure 9: Cumulative muon excesses from the direction of MknSOl (1° half- 
angle) . 



Cm 
O 

2 

I 

4 
3.5 



2.5 



1.5 



0.5 



ay 



CD 



01 



O) o 
O) o 
o 



MARKARIAN 421 



7jan93 14feb95 



v 



27 aug 97 



5 dec 98 




1000 



1500 



2000 



2500 



3000 



3500 



4000 



Julian Day (-2447500) 



Figure 10: Value of the quantity — log P (defined in the text) evaluated 
day- by-day from the direction of Mkn421 using a search cone with 1 "half- 
angle. The days with a value of — log P exceeding the probability value of 
10~^ that we have choosen as the attention level (indicated with the dashed 
line) are indicated by the arrows. 
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Figure 11: The same of the previous figure, but for the MknSOl source. 
No days exceeded the probability value of 10~^ that we have choosen as the 
attention level. 
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Figure 12: The cumulative distribution of the — log P quantity (see text) 
for all bins of the all-sky survey. 



